doseresNMA antidep/drnma.plot.1absolute.R

In htx-r/doseresNMA:

#!!!!! I used in max.dose- class instead of drug. Change when I plot it fro drug. Code it correctly
#!!!!! the rug is not added correctly with class - drug to class need to be changed maybe
#!!!!! add the rug properly or delete them from fig3(class)
# Plot: 
 # y1= absolute drug response, y2= placebo response VS
 # x=dose 
 # + rug of obseravtions
require(dplyr)
require(ggplot2)
# m=drnma_class
# data = antidep_class
# drug.lab=levels(antidep_class$class)
# nd=round(max(antidep_class$dose))
absolutePredplot <- function(m=m,
                             ref.lab='placebo',
                             data=antidep,
                             drug.lab=levels(antidep$drug),
                             nd=200,...){

  #-- Prepare data to plot
    # reference numeric index
    ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
    
    # 1. y-axis: p.drug 
    plotdata.drug <- m$BUGSoutput$summary %>% 
                 t() %>% 
                 data.frame() %>% 
                 select(starts_with('p.drug'))%>%
                 t()%>% 
                 data.frame() 
    # 2. add placebo data
    data.placebo <- m$BUGSoutput$summary%>% 
      t() %>% 
      data.frame() %>% 
      select(starts_with('p.placebo'))%>%
      t()%>% 
      data.frame() 
    plotdata <- plotdata.drug%>%
      mutate(X50.placebo=data.placebo$X50.,
             X2.5.placebo=data.placebo$X2.5.,
             X97.5.placebo=data.placebo$X97.5.)
    # 3. x-axis: doses per drug
     ndrugs <- length(drug.lab)-1 # number of drugs
     max.dose <- unlist(data%>%
      group_by(class)%>% #!!!!!!!!!!!!!! FIX class --> drug
      group_map(~round(max(.x$dose,na.rm = T))))[-ref.index]
     plotdata$dose <- unlist(
                            sapply(1:ndrugs, 
                                   function(i) seq(0,max.dose[i],l=nd),
                                   simplify = FALSE
                                   )
                            )
    
    
    # 3. labels: drug name
     plotdata$drug <- rep(drug.lab[drug.lab!=ref.lab],each=nd)
    
    # 4. rug: observed doses
     data0 <- data %>% # add a placebo indicator column
       group_by(studyid)%>%
       group_map(~mutate(.x,
                         is.placebo=ifelse(sum(.x$drug=='placebo')>0,0,NA)
       )
       ) # add 0 if we have a trial compare drug with placebo OR NA if not
     data0 <- data.frame(do.call(rbind,data0))
     
     plotdata$obs.dose <- unlist(sapply(1:ndrugs, 
                        function(i){
                          dose_plc_drug <- unique(unlist( # a vector of placebo(0/NA) and drug doses
                                                  data0[data0$drug%in%drug.lab[-ref.index][i],][c('is.placebo','dose')]
                                                    ))
                                                  
                          dose_plc_drug_nd <-rep(dose_plc_drug[!is.na(dose_plc_drug)],
                                                 l=nd) # repeat the vector until nd length to fit the data 
                          dose_plc_drug_nd
                                    }
                        ,simplify = F
                        )
                        )
  
  #-- To plot     
  # y axis limits
  ymax <- round(max(plotdata$X97.5.),1)
  ymin <- round(min(plotdata$X2.5.),1)
  
  # theme
  theme_set(
    theme_minimal() +
      theme(legend.position = "right",
            panel.background = element_rect(fill = 'snow1',colour = 'white'),
            axis.text.x = element_text(size=10),axis.text.y = element_text(size=10),
            axis.title.x=element_text(size=16,face = "bold"),axis.title.y=element_text(size=16,face = "bold"),
            strip.background =element_rect(fill="snow3"),
            strip.text.x = element_text(size = 14))
  )
  
  # plot
  g_X50 <- # the median line of drug and placebo
       ggplot(plotdata,aes(x=dose)) +       
       geom_line(aes(y=X50.,color='treatment'))+       # p.drug median
       geom_line(aes(y=X50.placebo,color='placebo'))+  # p.placebo: median
       scale_color_manual(values=c('steelblue','coral4'),name="")+ 
       guides(color = guide_legend(override.aes = list(size = 1.5)))
  
  g_crI <-  # add 95% CrI of p.drug and p.placebo
       g_X50+geom_smooth(                                    
       aes(x=dose,y=X50.,ymin=X2.5.,ymax=X97.5.),
       color='coral4',fill='coral1',
       data=plotdata, stat="identity")+ 
       geom_smooth(                                   
       aes(x=dose, y=X50.placebo, ymin=X2.5.placebo,ymax=X97.5.placebo),
       color='steelblue',fill='steelblue2',
       data=plotdata, stat="identity")+
       coord_cartesian(ylim = c(ymin, ymax))#+
  g_rug <- # add rug of the observed doses
           g_crI+geom_rug(#data=plotdata,
                   mapping=aes(x=obs.dose),
                   inherit.aes = F,
                   color='red',
                   length = unit(0.05, "npc"))
  
  g_split <- g_rug + facet_wrap(~drug,scales = 'free_x')
  
  g_labs <- g_split + ggplot2::labs(y="Predicted absolute response", x="Dose")
  
  g <- g_labs + ggplot2::scale_linetype_manual(name="")
  
  g
}